require(phyloseq)
Loading required package: phyloseq
require(tidyverse)
Loading required package: tidyverse
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5 ✓ purrr 0.3.4
✓ tibble 3.1.6 ✓ dplyr 1.0.8
✓ tidyr 1.2.0 ✓ stringr 1.4.0
✓ readr 2.1.2 ✓ forcats 0.5.1
Warning: package ‘tidyr’ was built under R version 4.1.2
Warning: package ‘readr’ was built under R version 4.1.2
Warning: package ‘dplyr’ was built under R version 4.1.2
── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
require(reshape2)
Loading required package: reshape2
Attaching package: ‘reshape2’
The following object is masked from ‘package:tidyr’:
smiths
require(dplyr)
require(ggplot2)
require(vegan)
Loading required package: vegan
Warning: package ‘vegan’ was built under R version 4.1.2
Loading required package: permute
Warning: package ‘permute’ was built under R version 4.1.2
Loading required package: lattice
This is vegan 2.6-2
Load data order, factors, and create a mode (chemical, hand, non-treated) column.
ps_dmn <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/DMN_ests_16S.Rdata")
sample_data(ps_dmn)$Herbicide <- factor(sample_data(ps_dmn)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_dmn)$herb_time<-paste(sample_data(ps_dmn)$Herbicide, sample_data(ps_dmn)$Time, sep = "_")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_dmn)$Mode<-sample_data(ps_dmn)$Herbicide
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")
sample_data(ps_dmn)$Mode<- as.factor(values[match(sample_data(ps_dmn)$Mode, index)])
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
sample_data(ps_dmn)$Herbicide <- as.factor(values[match(sample_data(ps_dmn)$Herbicide, index)])
ps_rare <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_rare_16S.Rdata")
sample_data(ps_rare)$Herbicide <- factor(sample_data(ps_rare)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_rare)$herb_time<-paste(sample_data(ps_rare)$Herbicide, sample_data(ps_rare)$Time, sep = "_")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")
sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
sample_data(ps_rare)$Herbicide <- as.factor(values[match(sample_data(ps_rare)$Herbicide, index)])
ps_trans <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_hel_trans_16S.Rdata")
sample_data(ps_trans)$Herbicide <- factor(sample_data(ps_trans)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_trans)$herb_time<-paste(sample_data(ps_trans)$Herbicide, sample_data(ps_trans)$Time, sep = "_")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_trans)$Mode<-sample_data(ps_trans)$Herbicide
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")
sample_data(ps_trans)$Mode<- as.factor(values[match(sample_data(ps_trans)$Mode, index)])
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
sample_data(ps_trans)$Herbicide <- as.factor(values[match(sample_data(ps_trans)$Herbicide, index)])
create alpha diversity tables
alpha_div_md$Herbicide
[1] Dicamba Dicamba Dicamba Handweeded Handweeded Handweeded Non-Treated Non-Treated
[9] Non-Treated Atrazine-Mesotrione Atrazine-Mesotrione Glyphosate Glyphosate Glyphosate Glyphosate Glyphosate
[17] Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded Handweeded Dicamba Dicamba Dicamba
[25] Non-Treated Non-Treated Non-Treated Non-Treated Non-Treated Non-Treated Glyphosate Glyphosate
[33] Glyphosate Dicamba Dicamba Dicamba Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded
[41] Handweeded Handweeded Glyphosate Glyphosate Glyphosate Handweeded Handweeded Non-Treated
[49] Non-Treated Non-Treated Dicamba Dicamba Dicamba Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione
[57] Dicamba Dicamba Handweeded Handweeded Handweeded Non-Treated Non-Treated Non-Treated
[65] Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Glyphosate Glyphosate Glyphosate Glyphosate Glyphosate
[73] Glyphosate Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded Handweeded Dicamba Dicamba
[81] Dicamba Non-Treated Non-Treated Non-Treated Non-Treated Non-Treated Glyphosate Glyphosate
[89] Dicamba Dicamba Dicamba Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded Handweeded
[97] Handweeded Glyphosate Glyphosate Handweeded Handweeded Handweeded Non-Treated Non-Treated
[105] Non-Treated Dicamba Dicamba Dicamba Atrazine-Mesotrione Atrazine-Mesotrione Dicamba Dicamba
[113] Dicamba Handweeded Handweeded Handweeded Non-Treated Non-Treated Atrazine-Mesotrione Atrazine-Mesotrione
[121] Atrazine-Mesotrione Glyphosate Glyphosate Glyphosate Glyphosate Glyphosate Atrazine-Mesotrione Atrazine-Mesotrione
[129] Atrazine-Mesotrione Handweeded Handweeded Handweeded Dicamba Dicamba Dicamba Non-Treated
[137] Non-Treated Non-Treated Non-Treated Non-Treated Glyphosate Glyphosate Glyphosate Dicamba
[145] Dicamba Dicamba Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded Handweeded Glyphosate
[153] Glyphosate Glyphosate Handweeded Handweeded Handweeded Non-Treated Non-Treated Non-Treated
[161] Dicamba Dicamba Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione
Levels: Non-Treated Handweeded Atrazine-Mesotrione Dicamba Glyphosate
Shannon Div plots - no significant differences among herbicide treatments at any of the three time points
ggplot(data = alpha_div_md, aes(Herbicide, Shannon, color= Herbicide)) + facet_grid(. ~ Time) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1) )
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_Shannon.pdf")
Saving 7.29 x 4.51 in image
aov_t1<-aov(Chao1 ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T1",])
plot(aov_t1$residuals)
summary(aov_t1)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 66537 16634 0.099 0.982
Residuals 51 8553352 167713
aov_t2<-aov(Chao1 ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T2",])
plot(aov_t2$residuals)
summary(aov_t2)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 652901 163225 0.779 0.545
Residuals 49 10272922 209651
aov_t3<-aov(Chao1 ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3",])
plot(aov_t3$residuals)
summary(aov_t3)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 433435 108359 0.465 0.761
Residuals 50 11641251 232825
remove outliers and rare reads with less than 2 total reads
ps_dmn <- subset_samples(ps_dmn, sample_names(ps_rare) != "G166SG")
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
ps_rare_sub<-prune_taxa(taxa_sums(ps_rare) > 2, ps_rare)
ps_trans_sub<-prune_taxa(taxa_sums(ps_trans) > 0.05, ps_trans)
ordinations and adonis testing with three separate objects (i.e., dmn, rarefied, transformed). Rare taxa are removed from rarefied and transfomred to sucessfully ordinate. At this point, the transformed data will not ordinate.
ord_dmn<-ordinate(physeq = ps_dmn, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit=1000)
Run 0 stress 0.1117941
Run 1 stress 0.1115908
... New best solution
... Procrustes: rmse 0.01112287 max resid 0.137285
Run 2 stress 0.1118138
... Procrustes: rmse 0.01224305 max resid 0.1486074
Run 3 stress 0.1117997
... Procrustes: rmse 0.01238183 max resid 0.1520467
Run 4 stress 0.1118133
... Procrustes: rmse 0.01165924 max resid 0.1445305
Run 5 stress 0.1117717
... Procrustes: rmse 0.004992797 max resid 0.05385656
Run 6 stress 0.1115916
... Procrustes: rmse 0.0003952995 max resid 0.003519588
... Similar to previous best
Run 7 stress 0.1117649
... Procrustes: rmse 0.004519199 max resid 0.05419881
Run 8 stress 0.1115917
... Procrustes: rmse 0.001064847 max resid 0.01208441
Run 9 stress 0.1117663
... Procrustes: rmse 0.00470262 max resid 0.0554819
Run 10 stress 0.1159638
Run 11 stress 0.1117987
... Procrustes: rmse 0.01229483 max resid 0.1505899
Run 12 stress 0.1118121
... Procrustes: rmse 0.01145765 max resid 0.1419449
Run 13 stress 0.1118606
... Procrustes: rmse 0.008173993 max resid 0.09623366
Run 14 stress 0.111798
... Procrustes: rmse 0.01121996 max resid 0.14002
Run 15 stress 0.1118121
... Procrustes: rmse 0.01340603 max resid 0.1643183
Run 16 stress 0.1115947
... Procrustes: rmse 0.001573293 max resid 0.01899481
Run 17 stress 0.1118639
... Procrustes: rmse 0.008635792 max resid 0.1014058
Run 18 stress 0.1115911
... Procrustes: rmse 0.0003000347 max resid 0.002739309
... Similar to previous best
Run 19 stress 0.1115954
... Procrustes: rmse 0.001243378 max resid 0.01498226
Run 20 stress 0.1117972
... Procrustes: rmse 0.01216683 max resid 0.1493349
*** Solution reached
ord_rare<-ordinate(physeq = ps_rare_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit=1000)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1959743
Run 1 stress 0.1955267
... New best solution
... Procrustes: rmse 0.01800803 max resid 0.1123914
Run 2 stress 0.1955394
... Procrustes: rmse 0.007006262 max resid 0.05765258
Run 3 stress 0.1958487
... Procrustes: rmse 0.01579267 max resid 0.1056545
Run 4 stress 0.1955945
... Procrustes: rmse 0.01295785 max resid 0.1494722
Run 5 stress 0.1955967
... Procrustes: rmse 0.01293347 max resid 0.1501704
Run 6 stress 0.1969796
Run 7 stress 0.1955973
... Procrustes: rmse 0.01337124 max resid 0.1534108
Run 8 stress 0.1955967
... Procrustes: rmse 0.01339074 max resid 0.1533846
Run 9 stress 0.1959732
... Procrustes: rmse 0.01729624 max resid 0.1080525
Run 10 stress 0.1956108
... Procrustes: rmse 0.01142375 max resid 0.1048192
Run 11 stress 0.1956435
... Procrustes: rmse 0.01260475 max resid 0.1537692
Run 12 stress 0.1955174
... New best solution
... Procrustes: rmse 0.008122313 max resid 0.05712863
Run 13 stress 0.1958337
... Procrustes: rmse 0.01090413 max resid 0.1036971
Run 14 stress 0.1967865
Run 15 stress 0.1955939
... Procrustes: rmse 0.01324124 max resid 0.1490696
Run 16 stress 0.196786
Run 17 stress 0.1955234
... Procrustes: rmse 0.007498469 max resid 0.05695487
Run 18 stress 0.195582
... Procrustes: rmse 0.01173812 max resid 0.07111958
Run 19 stress 0.1958207
... Procrustes: rmse 0.009376553 max resid 0.1035924
Run 20 stress 0.1967535
Run 21 stress 0.1956211
... Procrustes: rmse 0.01156007 max resid 0.1369555
Run 22 stress 0.195582
... Procrustes: rmse 0.01170171 max resid 0.07088219
Run 23 stress 0.1955424
... Procrustes: rmse 0.009362655 max resid 0.06422712
Run 24 stress 0.1954828
... New best solution
... Procrustes: rmse 0.006018466 max resid 0.04375112
Run 25 stress 0.1955168
... Procrustes: rmse 0.00608034 max resid 0.04463714
Run 26 stress 0.1965865
Run 27 stress 0.196778
Run 28 stress 0.1955925
... Procrustes: rmse 0.01221595 max resid 0.1482877
Run 29 stress 0.1955944
... Procrustes: rmse 0.01222552 max resid 0.1482843
Run 30 stress 0.1958206
... Procrustes: rmse 0.01151099 max resid 0.1048252
Run 31 stress 0.1956224
... Procrustes: rmse 0.01254164 max resid 0.134756
Run 32 stress 0.1971977
Run 33 stress 0.1969766
Run 34 stress 0.195833
... Procrustes: rmse 0.01382711 max resid 0.1051923
Run 35 stress 0.1954885
... Procrustes: rmse 0.001224376 max resid 0.01288981
Run 36 stress 0.1954874
... Procrustes: rmse 0.001327997 max resid 0.01348466
Run 37 stress 0.1967515
Run 38 stress 0.1958228
... Procrustes: rmse 0.01137911 max resid 0.1048475
Run 39 stress 0.197019
Run 40 stress 0.1958336
... Procrustes: rmse 0.01400106 max resid 0.1052911
Run 41 stress 0.1955898
... Procrustes: rmse 0.01087936 max resid 0.1192967
Run 42 stress 0.1954861
... Procrustes: rmse 0.001352308 max resid 0.01107595
Run 43 stress 0.1984585
Run 44 stress 0.1955182
... Procrustes: rmse 0.005903524 max resid 0.04272373
Run 45 stress 0.1970057
Run 46 stress 0.1965857
Run 47 stress 0.1974123
Run 48 stress 0.1956422
... Procrustes: rmse 0.01285919 max resid 0.1537516
Run 49 stress 0.1956489
... Procrustes: rmse 0.01316701 max resid 0.1561775
Run 50 stress 0.195849
... Procrustes: rmse 0.01454653 max resid 0.1054922
Run 51 stress 0.195821
... Procrustes: rmse 0.01228135 max resid 0.105193
Run 52 stress 0.1978668
Run 53 stress 0.1958338
... Procrustes: rmse 0.01282964 max resid 0.10528
Run 54 stress 0.1955268
... Procrustes: rmse 0.003679333 max resid 0.03000048
Run 55 stress 0.1971375
Run 56 stress 0.1956175
... Procrustes: rmse 0.01127198 max resid 0.1151186
Run 57 stress 0.1956146
... Procrustes: rmse 0.01061111 max resid 0.1067233
Run 58 stress 0.1958342
... Procrustes: rmse 0.01302074 max resid 0.1053564
Run 59 stress 0.1972164
Run 60 stress 0.196753
Run 61 stress 0.197357
Run 62 stress 0.1956426
... Procrustes: rmse 0.01282273 max resid 0.1539347
Run 63 stress 0.1956151
... Procrustes: rmse 0.01074056 max resid 0.1084906
Run 64 stress 0.195649
... Procrustes: rmse 0.01317392 max resid 0.1561167
Run 65 stress 0.1959754
... Procrustes: rmse 0.01604176 max resid 0.1104391
Run 66 stress 0.1954887
... Procrustes: rmse 0.001712985 max resid 0.01764067
Run 67 stress 0.1969292
Run 68 stress 0.1967917
Run 69 stress 0.1955019
... Procrustes: rmse 0.003005991 max resid 0.01927756
Run 70 stress 0.196758
Run 71 stress 0.1972331
Run 72 stress 0.1971941
Run 73 stress 0.1954839
... Procrustes: rmse 0.000538266 max resid 0.005283253
... Similar to previous best
*** Solution reached
full_ord_rare<-ggordiplots::gg_ordiplot(ord = ord_rare, groups = data.frame(sample_data(ps_rare))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
full_ord_rare$plot + theme_classic()
ord_transformed<-ordinate(physeq = ps_trans_sub, method = "NMDS", distance = "bray", trymax= 300, maxit=1000)
Run 0 stress 0.07436383
Run 1 stress 0.07411606
... New best solution
... Procrustes: rmse 0.03117035 max resid 0.2498244
Run 2 stress 0.07773971
Run 3 stress 0.07217597
... New best solution
... Procrustes: rmse 0.05313851 max resid 0.4765044
Run 4 stress 0.07455535
Run 5 stress 0.07413861
Run 6 stress 0.07821068
Run 7 stress 0.07475795
Run 8 stress 0.07634371
Run 9 stress 0.07390696
Run 10 stress 0.07475868
Run 11 stress 0.07371342
Run 12 stress 0.07653065
Run 13 stress 0.07356773
Run 14 stress 0.07214992
... New best solution
... Procrustes: rmse 0.001443487 max resid 0.01837558
Run 15 stress 0.07478238
Run 16 stress 0.07401673
Run 17 stress 0.07564706
Run 18 stress 0.07376447
Run 19 stress 0.0763201
Run 20 stress 0.08031617
Run 21 stress 0.07350122
Run 22 stress 0.07472829
Run 23 stress 0.07455301
Run 24 stress 0.07224441
... Procrustes: rmse 0.002650218 max resid 0.02493022
Run 25 stress 0.07343235
Run 26 stress 0.08340946
Run 27 stress 0.07371382
Run 28 stress 0.07809551
Run 29 stress 0.07565057
Run 30 stress 0.07604995
Run 31 stress 0.07400801
Run 32 stress 0.07455409
Run 33 stress 0.08407341
Run 34 stress 0.08405253
Run 35 stress 0.08494465
Run 36 stress 0.07207254
... New best solution
... Procrustes: rmse 0.02923975 max resid 0.2677599
Run 37 stress 0.07433151
Run 38 stress 0.08031236
Run 39 stress 0.07401128
Run 40 stress 0.07630296
Run 41 stress 0.07504806
Run 42 stress 0.07353092
Run 43 stress 0.07210215
... Procrustes: rmse 0.001996828 max resid 0.0169429
Run 44 stress 0.07290392
Run 45 stress 0.07456734
Run 46 stress 0.4158745
Run 47 stress 0.07468176
Run 48 stress 0.07372598
Run 49 stress 0.07314001
Run 50 stress 0.07431572
Run 51 stress 0.07217868
... Procrustes: rmse 0.02929952 max resid 0.2671141
Run 52 stress 0.07401337
Run 53 stress 0.07419621
Run 54 stress 0.0761684
Run 55 stress 0.07390179
Run 56 stress 0.08495755
Run 57 stress 0.07457704
Run 58 stress 0.07372457
Run 59 stress 0.07411643
Run 60 stress 0.07810726
Run 61 stress 0.07435209
Run 62 stress 0.07565289
Run 63 stress 0.07399032
Run 64 stress 0.08507366
Run 65 stress 0.07419598
Run 66 stress 0.08221774
Run 67 stress 0.07388912
Run 68 stress 0.08249078
Run 69 stress 0.07401304
Run 70 stress 0.07388093
Run 71 stress 0.07399428
Run 72 stress 0.08033213
Run 73 stress 0.07371329
Run 74 stress 0.07379163
Run 75 stress 0.08424884
Run 76 stress 0.07382315
Run 77 stress 0.07495688
Run 78 stress 0.07370127
Run 79 stress 0.07475771
Run 80 stress 0.07215002
... Procrustes: rmse 0.02929661 max resid 0.2675789
Run 81 stress 0.07453389
Run 82 stress 0.07291297
Run 83 stress 0.08293632
Run 84 stress 0.07720927
Run 85 stress 0.07438212
Run 86 stress 0.07390787
Run 87 stress 0.07434478
Run 88 stress 0.07433927
Run 89 stress 0.08539586
Run 90 stress 0.08382492
Run 91 stress 0.07398725
Run 92 stress 0.07407376
Run 93 stress 0.07353074
Run 94 stress 0.07432943
Run 95 stress 0.07207262
... Procrustes: rmse 0.0004636339 max resid 0.004175306
... Similar to previous best
*** Solution reached
Adonis testing of herbicide treatments by time point
ps_adonis<-function(physeq){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
print(anova(vegan::betadisper(physeq_dist, md_tab$Herbicide)))
vegan::adonis(physeq_dist ~ Herbicide * Time + Total_Weed_Veg , data = md_tab, permutations = 1000)
}
remove one sample with no vegetation measurement.
ps_rare_sub_57<-subset_samples(ps_rare_sub, sample_names(ps_rare_sub) != "G065SG")
ps_adonis(ps_rare_sub_57)
ps_dmn_57<-subset_samples(ps_dmn, sample_names(ps_dmn) != "G065SG")
ps_adonis(ps_dmn_57)
Ordination plots DMN
ord_t1_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Run 0 stress 0.0860376
Run 1 stress 0.08603806
... Procrustes: rmse 0.0001239324 max resid 0.0007295483
... Similar to previous best
Run 2 stress 0.0983852
Run 3 stress 0.08674017
Run 4 stress 0.09974458
Run 5 stress 0.08662921
Run 6 stress 0.08662975
Run 7 stress 0.08674026
Run 8 stress 0.08603767
... Procrustes: rmse 2.019457e-05 max resid 0.0001271086
... Similar to previous best
Run 9 stress 0.1024486
Run 10 stress 0.0866297
Run 11 stress 0.0866588
Run 12 stress 0.0860374
... New best solution
... Procrustes: rmse 0.0008541763 max resid 0.004577351
... Similar to previous best
Run 13 stress 0.08674053
Run 14 stress 0.08662909
Run 15 stress 0.0997284
Run 16 stress 0.08603683
... New best solution
... Procrustes: rmse 0.0002639494 max resid 0.001218145
... Similar to previous best
Run 17 stress 0.08674004
Run 18 stress 0.08603779
... Procrustes: rmse 0.0006511547 max resid 0.003685441
... Similar to previous best
Run 19 stress 0.08673995
Run 20 stress 0.08673991
*** Solution reached
T1_dmn<-ggordiplots::gg_ordiplot(ord = ord_t1_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_dmn$plot + theme_classic() + xlim(-0.4, 0.4) + ylim(-0.4, 0.4)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T1.pdf")
Saving 7.29 x 4.51 in image
ord_t2_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Run 0 stress 0.1104159
Run 1 stress 0.1104153
... New best solution
... Procrustes: rmse 0.0003571464 max resid 0.002189988
... Similar to previous best
Run 2 stress 0.1095397
... New best solution
... Procrustes: rmse 0.03309861 max resid 0.2348476
Run 3 stress 0.1096257
... Procrustes: rmse 0.005714701 max resid 0.03065104
Run 4 stress 0.1095261
... New best solution
... Procrustes: rmse 0.0006951312 max resid 0.003638642
... Similar to previous best
Run 5 stress 0.1104161
Run 6 stress 0.1122104
Run 7 stress 0.1095262
... Procrustes: rmse 0.0001081361 max resid 0.0004776
... Similar to previous best
Run 8 stress 0.1122142
Run 9 stress 0.1095263
... Procrustes: rmse 0.0006652566 max resid 0.003388802
... Similar to previous best
Run 10 stress 0.1095263
... Procrustes: rmse 0.0005610351 max resid 0.003406831
... Similar to previous best
Run 11 stress 0.1108789
Run 12 stress 0.1095267
... Procrustes: rmse 0.0008601728 max resid 0.004694253
... Similar to previous best
Run 13 stress 0.1104149
Run 14 stress 0.1122102
Run 15 stress 0.1104149
Run 16 stress 0.1108808
Run 17 stress 0.1104147
Run 18 stress 0.1105647
Run 19 stress 0.1108803
Run 20 stress 0.1108777
*** Solution reached
T2_dmn<-ggordiplots::gg_ordiplot(ord = ord_t2_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_dmn$plot + theme_classic()+ xlim(-0.4, 0.4) + ylim(-0.4, 0.4)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T2.pdf")
Saving 7.29 x 4.51 in image
ord_t3_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Run 0 stress 0.09546215
Run 1 stress 0.09706147
Run 2 stress 0.09546159
... New best solution
... Procrustes: rmse 0.001277492 max resid 0.007259895
... Similar to previous best
Run 3 stress 0.0955498
... Procrustes: rmse 0.03810422 max resid 0.2424179
Run 4 stress 0.09550233
... Procrustes: rmse 0.03778166 max resid 0.2444873
Run 5 stress 0.09546332
... Procrustes: rmse 0.001100842 max resid 0.005583535
... Similar to previous best
Run 6 stress 0.09732995
Run 7 stress 0.09546077
... New best solution
... Procrustes: rmse 0.0007947735 max resid 0.00456194
... Similar to previous best
Run 8 stress 0.09546428
... Procrustes: rmse 0.001848437 max resid 0.009343808
... Similar to previous best
Run 9 stress 0.09712155
Run 10 stress 0.09549825
... Procrustes: rmse 0.03714564 max resid 0.2419844
Run 11 stress 0.09549948
... Procrustes: rmse 0.0373624 max resid 0.2433029
Run 12 stress 0.09546252
... Procrustes: rmse 0.0005053034 max resid 0.002683928
... Similar to previous best
Run 13 stress 0.09569417
... Procrustes: rmse 0.03993933 max resid 0.2524409
Run 14 stress 0.09546543
... Procrustes: rmse 0.00135589 max resid 0.007233676
... Similar to previous best
Run 15 stress 0.0971071
Run 16 stress 0.09549812
... Procrustes: rmse 0.03722254 max resid 0.2427672
Run 17 stress 0.0970629
Run 18 stress 0.09546061
... New best solution
... Procrustes: rmse 0.000334493 max resid 0.00147327
... Similar to previous best
Run 19 stress 0.09570193
... Procrustes: rmse 0.04006183 max resid 0.2508731
Run 20 stress 0.09549861
... Procrustes: rmse 0.03739865 max resid 0.2432734
*** Solution reached
T3_dmn<-ggordiplots::gg_ordiplot(ord = ord_t3_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_dmn$plot + theme_classic()+ xlim(-0.4, 0.4) + ylim(-0.4, 0.4)
Warning: Removed 1 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T3.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing missing values (geom_point).
Ordination plots rarefied
ord_t1_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1810491
Run 1 stress 0.1812643
... Procrustes: rmse 0.02583058 max resid 0.1092403
Run 2 stress 0.1814142
... Procrustes: rmse 0.04648293 max resid 0.1540837
Run 3 stress 0.1824085
Run 4 stress 0.1854315
Run 5 stress 0.1817696
Run 6 stress 0.1813795
... Procrustes: rmse 0.01903201 max resid 0.1229428
Run 7 stress 0.1873045
Run 8 stress 0.1810319
... New best solution
... Procrustes: rmse 0.03213373 max resid 0.1238225
Run 9 stress 0.1815254
... Procrustes: rmse 0.07330822 max resid 0.260521
Run 10 stress 0.181525
... Procrustes: rmse 0.07331803 max resid 0.2605261
Run 11 stress 0.1815251
... Procrustes: rmse 0.07331039 max resid 0.2605071
Run 12 stress 0.1812636
... Procrustes: rmse 0.02331833 max resid 0.1081706
Run 13 stress 0.1873031
Run 14 stress 0.1809659
... New best solution
... Procrustes: rmse 0.04119361 max resid 0.1474758
Run 15 stress 0.1822765
Run 16 stress 0.1864799
Run 17 stress 0.1815253
Run 18 stress 0.1837641
Run 19 stress 0.1814576
... Procrustes: rmse 0.05589707 max resid 0.2469422
Run 20 stress 0.1836188
Run 21 stress 0.1819126
Run 22 stress 0.1821857
Run 23 stress 0.1810495
... Procrustes: rmse 0.01875174 max resid 0.1152783
Run 24 stress 0.1822761
Run 25 stress 0.1808423
... New best solution
... Procrustes: rmse 0.03509173 max resid 0.1424772
Run 26 stress 0.1815222
Run 27 stress 0.1814211
Run 28 stress 0.1881278
Run 29 stress 0.1831937
Run 30 stress 0.1815254
Run 31 stress 0.1815254
Run 32 stress 0.1877224
Run 33 stress 0.1823821
Run 34 stress 0.1816035
Run 35 stress 0.1828841
Run 36 stress 0.1811594
... Procrustes: rmse 0.02182468 max resid 0.1078205
Run 37 stress 0.181438
Run 38 stress 0.182276
Run 39 stress 0.1824766
Run 40 stress 0.1875966
Run 41 stress 0.182182
Run 42 stress 0.1878771
Run 43 stress 0.1872815
Run 44 stress 0.1885431
Run 45 stress 0.1821561
Run 46 stress 0.1815255
Run 47 stress 0.1816037
Run 48 stress 0.1810502
... Procrustes: rmse 0.02488083 max resid 0.1084949
Run 49 stress 0.1842028
Run 50 stress 0.1831935
Run 51 stress 0.181416
Run 52 stress 0.1825128
Run 53 stress 0.1821252
Run 54 stress 0.1882
Run 55 stress 0.1858117
Run 56 stress 0.1910572
Run 57 stress 0.1814036
Run 58 stress 0.1815253
Run 59 stress 0.1837646
Run 60 stress 0.1855356
Run 61 stress 0.1886207
Run 62 stress 0.1811714
... Procrustes: rmse 0.02206901 max resid 0.1101655
Run 63 stress 0.1813283
... Procrustes: rmse 0.02669062 max resid 0.1245773
Run 64 stress 0.1815262
Run 65 stress 0.181948
Run 66 stress 0.1917236
Run 67 stress 0.1868432
Run 68 stress 0.1828928
Run 69 stress 0.1808185
... New best solution
... Procrustes: rmse 0.006158855 max resid 0.03531337
Run 70 stress 0.1819127
Run 71 stress 0.1815252
Run 72 stress 0.1886606
Run 73 stress 0.186167
Run 74 stress 0.180948
... Procrustes: rmse 0.009886941 max resid 0.05098863
Run 75 stress 0.1816035
Run 76 stress 0.1935237
Run 77 stress 0.1812551
... Procrustes: rmse 0.01336276 max resid 0.08156214
Run 78 stress 0.1822264
Run 79 stress 0.1814582
Run 80 stress 0.1814144
Run 81 stress 0.1828653
Run 82 stress 0.1823118
Run 83 stress 0.1978016
Run 84 stress 0.1810713
... Procrustes: rmse 0.0369974 max resid 0.1455491
Run 85 stress 0.1808181
... New best solution
... Procrustes: rmse 0.0002833112 max resid 0.0008111158
... Similar to previous best
*** Solution reached
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare_plot<-T1_rare$plot + theme_classic() + xlim(-0.4, 0.4) + ylim(-0.4, 0.4) + guides(color=guide_legend("Treatment")) + xlab("NMDS 1") + ylab("NMDS 2")
T1_rare_plot
library(cowplot)
my_legend <- get_legend(T1_rare_plot)
library(ggpubr)
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordinationlegend.pdf")
Saving 7.29 x 4.51 in image
T1_rare_plot<-T1_rare_plot + theme(legend.position = "none")
T1_rare_plot
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T1.pdf")
Saving 7.29 x 4.51 in image
ord_t2_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1534021
Run 1 stress 0.1495833
... New best solution
... Procrustes: rmse 0.06591139 max resid 0.2432552
Run 2 stress 0.1529359
Run 3 stress 0.149299
... New best solution
... Procrustes: rmse 0.01311101 max resid 0.05357392
Run 4 stress 0.1498585
Run 5 stress 0.153965
Run 6 stress 0.1535223
Run 7 stress 0.1501546
Run 8 stress 0.1508771
Run 9 stress 0.1500086
Run 10 stress 0.154203
Run 11 stress 0.1496264
... Procrustes: rmse 0.04805853 max resid 0.2154931
Run 12 stress 0.1546731
Run 13 stress 0.1496216
... Procrustes: rmse 0.04712721 max resid 0.2093185
Run 14 stress 0.15183
Run 15 stress 0.1496336
... Procrustes: rmse 0.03080855 max resid 0.1376397
Run 16 stress 0.1496321
... Procrustes: rmse 0.03047098 max resid 0.13453
Run 17 stress 0.1500846
Run 18 stress 0.1495175
... Procrustes: rmse 0.02694526 max resid 0.1366575
Run 19 stress 0.1517768
Run 20 stress 0.1498159
Run 21 stress 0.1525828
Run 22 stress 0.1496191
... Procrustes: rmse 0.04687338 max resid 0.205771
Run 23 stress 0.1558324
Run 24 stress 0.1558559
Run 25 stress 0.1517661
Run 26 stress 0.1493063
... Procrustes: rmse 0.00273002 max resid 0.01544446
Run 27 stress 0.1543836
Run 28 stress 0.1493533
... Procrustes: rmse 0.02753152 max resid 0.1459909
Run 29 stress 0.1497401
... Procrustes: rmse 0.02124904 max resid 0.09580409
Run 30 stress 0.1536117
Run 31 stress 0.152943
Run 32 stress 0.153367
Run 33 stress 0.1502911
Run 34 stress 0.1492996
... Procrustes: rmse 0.009094635 max resid 0.05431342
Run 35 stress 0.1500939
Run 36 stress 0.157447
Run 37 stress 0.1501923
Run 38 stress 0.1503452
Run 39 stress 0.1499623
Run 40 stress 0.1565076
Run 41 stress 0.1552341
Run 42 stress 0.1493037
... Procrustes: rmse 0.01222299 max resid 0.07016078
Run 43 stress 0.1542397
Run 44 stress 0.1504335
Run 45 stress 0.1502563
Run 46 stress 0.1493242
... Procrustes: rmse 0.006378492 max resid 0.03651169
Run 47 stress 0.1493052
... Procrustes: rmse 0.002363587 max resid 0.01335982
Run 48 stress 0.1496238
... Procrustes: rmse 0.04801997 max resid 0.2142248
Run 49 stress 0.1493332
... Procrustes: rmse 0.007566515 max resid 0.04258916
Run 50 stress 0.1508837
Run 51 stress 0.1522037
Run 52 stress 0.1492999
... Procrustes: rmse 0.01056149 max resid 0.06267914
Run 53 stress 0.1504227
Run 54 stress 0.1495976
... Procrustes: rmse 0.01568511 max resid 0.07203582
Run 55 stress 0.154981
Run 56 stress 0.1541634
Run 57 stress 0.1492997
... Procrustes: rmse 0.001503889 max resid 0.006022071
... Similar to previous best
*** Solution reached
T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_rare_plot<-T2_rare$plot + theme_classic() + xlim(-0.4, 0.4) + ylim(-0.4, 0.4)+ theme(legend.position = "none") + xlab("NMDS 1") + ylab("NMDS 2")
T2_rare_plot
Warning: Removed 3 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T2.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 3 rows containing missing values (geom_point).
#G166SG identified as outlier based on plots with it included. Removed to create plot.
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
ord_t3_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1676017
Run 1 stress 0.1674408
... New best solution
... Procrustes: rmse 0.04160458 max resid 0.2358723
Run 2 stress 0.1676017
... Procrustes: rmse 0.04176362 max resid 0.2370166
Run 3 stress 0.1676781
... Procrustes: rmse 0.03928737 max resid 0.234611
Run 4 stress 0.1674409
... Procrustes: rmse 0.0002702533 max resid 0.00115854
... Similar to previous best
Run 5 stress 0.1674412
... Procrustes: rmse 0.0009690672 max resid 0.005153489
... Similar to previous best
Run 6 stress 0.1674405
... New best solution
... Procrustes: rmse 0.0005420413 max resid 0.002292957
... Similar to previous best
Run 7 stress 0.1674412
... Procrustes: rmse 0.00046538 max resid 0.002452812
... Similar to previous best
Run 8 stress 0.1752381
Run 9 stress 0.1674408
... Procrustes: rmse 0.0002149259 max resid 0.0009482109
... Similar to previous best
Run 10 stress 0.1674412
... Procrustes: rmse 0.0005201256 max resid 0.002801695
... Similar to previous best
Run 11 stress 0.1674417
... Procrustes: rmse 0.0006125434 max resid 0.002951462
... Similar to previous best
Run 12 stress 0.167681
... Procrustes: rmse 0.03993784 max resid 0.2358893
Run 13 stress 0.1674407
... Procrustes: rmse 0.0003510634 max resid 0.001544427
... Similar to previous best
Run 14 stress 0.1676786
... Procrustes: rmse 0.03985713 max resid 0.235582
Run 15 stress 0.1674409
... Procrustes: rmse 0.0005786015 max resid 0.003398217
... Similar to previous best
Run 16 stress 0.1676795
... Procrustes: rmse 0.03999965 max resid 0.2356988
Run 17 stress 0.1674408
... Procrustes: rmse 0.000393707 max resid 0.001767988
... Similar to previous best
Run 18 stress 0.167441
... Procrustes: rmse 0.0003750548 max resid 0.001801765
... Similar to previous best
Run 19 stress 0.1674412
... Procrustes: rmse 0.0005055888 max resid 0.002397934
... Similar to previous best
Run 20 stress 0.1676016
... Procrustes: rmse 0.0417145 max resid 0.2369253
*** Solution reached
T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_rare_plot<-T3_rare$plot + theme_classic() + xlim(-0.4, 0.4) + ylim(-0.4, 0.4)+ theme(legend.position = "none") + xlab("NMDS 1") + ylab("NMDS 2")
T3_rare_plot
Warning: Removed 1 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T3.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing missing values (geom_point).
library(ggpubr)
ggarrange(T1_rare_plot, T2_rare_plot, T3_rare_plot, ncol = 1)
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 1 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta.pdf", width = 10, height = 10)
CAP ordination plots rarefied - not used
t1_dist <- distance(subset_samples(ps_rare, Time=="T1"), method="bray") #get wUnifrac and save
t1_table<-as.matrix(dist(t1_dist)) #transform wUnifrac index
ord_t1_rare_cap <- capscale(t1_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T1"))))
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T1_cap.pdf")
t2_dist <- distance(subset_samples(ps_rare, Time=="T2"), method="bray") #get wUnifrac and save
t2_table<-as.matrix(dist(t2_dist)) #transform wUnifrac index
ord_t2_rare_cap <- capscale(t2_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T2"))))
T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T2_cap.pdf")
#G166SG identified as outlier based on plots with it included. Removed to create plot.
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
t3_dist <- distance(subset_samples(ps_rare, Time=="T3"), method="bray") #get wUnifrac and save
t3_table<-as.matrix(dist(t3_dist)) #transform wUnifrac index
ord_t3_rare_cap <- capscale(t3_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T3"))))
T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T3_cap.pdf")
ggplot_build(T3_rare$plot)$data
Pairwise adonis testing
ps_pairwiseadonis<-function(physeq){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
pairwiseAdonis::pairwise.adonis(x = physeq_dist, factors = md_tab$Herbicide, p.adjust.m = "none", perm = 1000)
}
ps_t1<-subset_samples(ps_rare_sub, Time == "T1")
ps_t1<-prune_taxa(taxa_sums(ps_t1) > 2, ps_t1)
ps_t2<-subset_samples(ps_rare_sub, Time == "T2")
ps_t2<-prune_taxa(taxa_sums(ps_t2) > 2, ps_t2)
ps_t3<-subset_samples(ps_rare_sub, Time == "T3")
ps_t3<-prune_taxa(taxa_sums(ps_t3) > 2, ps_t3)
ps_pairwiseadonis(ps_t1)
ps_pairwiseadonis(ps_t2)
ps_pairwiseadonis(ps_t3)
Pairwise betadispr by treatment, time and mode
ps_betadispr<-function(physeq, groupingvar = "Groupingvar"){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
mod<-vegan::betadisper(physeq_dist, md_tab[,groupingvar])
## Perform test
print(anova(mod))
## Permutation test for F
pmod <- vegan::permutest(mod, permutations = 1000, pairwise = TRUE)
print(pmod)
print(boxplot(mod))
}
permute test of disperson
ps_betadispr(subset_samples(ps_rare_sub, Time == "T1"), groupingvar = "Mode")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T2"), groupingvar = "Mode")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T3"), groupingvar = "Mode")
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Chemical"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Non-Treated"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Hand"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T1"), groupingvar = "Herbicide")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T2"), groupingvar = "Herbicide")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T3"), groupingvar = "Herbicide")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Glyphosate"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Atrazine-Mesotrione"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Dicamba"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Handweeded"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Non-Treated"), groupingvar = "Time")
ps_betadispr(ps_rare_sub, groupingvar = "Herbicide")
ps_betadispr(ps_rare_sub, groupingvar = "Mode")
ps_betadispr(ps_rare_sub, groupingvar = "Time")
box and whisker plots of pairwise distance within group distances
#remotes::install_github("antonioggsousa/micrUBIfuns")
library(micrUBIfuns)
T1_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Herbicide")
T1_beta_plot <- T1_beta$plot
T1_beta_plot <- T1_beta_plot + theme_classic()+ guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75)
T1_beta_plot
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
my_legend <- get_legend(T1_beta_plot)
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_beta_legend.pdf")
Saving 7.29 x 4.51 in image
T1_beta_plot<-T1_beta_plot+ theme(legend.position = "none")
T1_beta_plot
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
T1_beta_df<- T1_beta$data
T1_betamod<-aov(formula = beta_div_value ~ group ,data = T1_beta_df)
summary(T1_betamod)
Df Sum Sq Mean Sq F value Pr(>F)
group 4 0.02391 0.005978 5.4 0.000334 ***
Residuals 282 0.31223 0.001107
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(x = T1_betamod, which = "group")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ group, data = T1_beta_df)
$group
diff lwr upr p adj
Dicamba-Atrazine-Mesotrione -0.001507071 -0.018186144 0.015172002 0.9991606
Glyphosate-Atrazine-Mesotrione 0.001896970 -0.015523754 0.019317693 0.9982525
Handweeded-Atrazine-Mesotrione -0.023576431 -0.041939486 -0.005213376 0.0044684
Non-Treated-Atrazine-Mesotrione -0.013193939 -0.029873012 0.003485134 0.1934686
Glyphosate-Dicamba 0.003404040 -0.013275033 0.020083113 0.9805780
Handweeded-Dicamba -0.022069360 -0.039730381 -0.004408340 0.0061800
Non-Treated-Dicamba -0.011686869 -0.027589741 0.004216003 0.2601887
Handweeded-Glyphosate -0.025473401 -0.043836456 -0.007110346 0.0015994
Non-Treated-Glyphosate -0.015090909 -0.031769982 0.001588164 0.0971493
Non-Treated-Handweeded 0.010382492 -0.007278529 0.028043512 0.4896459
T2_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Herbicide")
T2_beta_plot <- T2_beta$plot
T2_beta_plot <- T2_beta_plot+ theme_classic() + theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75)
T2_beta_plot
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
T2_beta_df<- T2_beta$data
T2_betamod<-aov(formula = beta_div_value ~ group ,data = T2_beta_df)
summary(T2_betamod)
Df Sum Sq Mean Sq F value Pr(>F)
group 4 0.0328 0.008212 5.471 0.000303 ***
Residuals 260 0.3902 0.001501
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(x = T2_betamod, which = "group")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ group, data = T2_beta_df)
$group
diff lwr upr p adj
Dicamba-Atrazine-Mesotrione 0.004587879 -0.015705828 0.024881585 0.9716343
Glyphosate-Atrazine-Mesotrione -0.011421549 -0.032812994 0.009969896 0.5850476
Handweeded-Atrazine-Mesotrione 0.007709091 -0.012584615 0.028002797 0.8348219
Non-Treated-Atrazine-Mesotrione -0.022036364 -0.042330070 -0.001742657 0.0257655
Glyphosate-Dicamba -0.016009428 -0.037400872 0.005382017 0.2426725
Handweeded-Dicamba 0.003121212 -0.017172494 0.023414919 0.9933115
Non-Treated-Dicamba -0.026624242 -0.046917949 -0.006330536 0.0034255
Handweeded-Glyphosate 0.019130640 -0.002260805 0.040522085 0.1039333
Non-Treated-Glyphosate -0.010614815 -0.032006260 0.010776630 0.6518171
Non-Treated-Handweeded -0.029745455 -0.050039161 -0.009451748 0.0007049
T3_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Herbicide")
T3_beta$plot #+ scale_color_manual(values = c("#F8766D", "#A3A500", "#00BF7D", "#00B0F6", "#E76BF3")) +
T3_beta_plot <- T3_beta$plot
T3_beta_plot <- T3_beta_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("")
T3_beta_plot <-T3_beta_plot + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
T3_beta_df<- T3_beta$data
T3_betamod<-aov(formula = beta_div_value ~ group ,data = T3_beta_df)
summary(T3_betamod)
Df Sum Sq Mean Sq F value Pr(>F)
group 4 0.0477 0.011924 10.39 7.93e-08 ***
Residuals 261 0.2994 0.001147
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(x = T3_betamod, which = "group")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ group, data = T3_beta_df)
$group
diff lwr upr p adj
Dicamba-Atrazine-Mesotrione 0.0183969697 0.001410694 0.035383246 0.0263568
Glyphosate-Atrazine-Mesotrione -0.0007666667 -0.018752976 0.017219643 0.9999574
Handweeded-Atrazine-Mesotrione 0.0276454545 0.010659179 0.044631731 0.0001130
Non-Treated-Atrazine-Mesotrione 0.0315000000 0.013513690 0.049486310 0.0000250
Glyphosate-Dicamba -0.0191636364 -0.037864911 -0.000462362 0.0415656
Handweeded-Dicamba 0.0092484848 -0.008493102 0.026990072 0.6075977
Non-Treated-Dicamba 0.0131030303 -0.005598244 0.031804305 0.3068486
Handweeded-Glyphosate 0.0284121212 0.009710847 0.047113396 0.0003917
Non-Treated-Glyphosate 0.0322666667 0.012652605 0.051880729 0.0000917
Non-Treated-Handweeded 0.0038545455 -0.014846729 0.022555820 0.9798083
library(ggpubr)
ggarrange(T1_beta_plot, T2_beta_plot, T3_beta_plot, ncol = 1)
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta.pdf", width = 5, height = 10)
Examination of dissimliarity across time points by treatment and then again by all chemical treatments combined.
T1_beta_df$Time<-"T1"
T2_beta_df$Time<-"T2"
T3_beta_df$Time<-"T3"
beta_div_T1_T2_T3 <- rbind(T1_beta_df, T2_beta_df, T3_beta_df)
beta_anova<-function(data, Herbicide = "Herbicide"){
df_sub<- data %>% filter(group == Herbicide)
mod<-aov(beta_div_value ~ Time, data = df_sub)
print(summary(mod))
print(TukeyHSD(mod, "Time"))
boxplot(df_sub$beta_div_value ~ df_sub$Time)
}
beta_anova(beta_div_T1_T2_T3, Herbicide = "Non-Treated")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.02547 0.01274 17.94 9.08e-08 ***
Residuals 163 0.11573 0.00071
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 -0.01621212 -0.027718988 -0.004705255 0.0030431
T3-T1 0.01578788 0.003603568 0.027972190 0.0071575
T3-T2 0.03200000 0.019331357 0.044668643 0.0000000
beta_anova(beta_div_T1_T2_T3, Herbicide = "Handweeded")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.01713 0.008567 4.195 0.0168 *
Residuals 152 0.31041 0.002042
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 0.02391582 0.0024158024 0.04541585 0.0251783
T3-T1 0.02231582 0.0008158024 0.04381585 0.0399463
T3-T2 -0.00160000 -0.0219967123 0.01879671 0.9811772
beta_anova(beta_div_T1_T2_T3, Herbicide = "Dicamba")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.00273 0.001366 0.977 0.378
Residuals 173 0.24171 0.001397
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 -0.001274747 -0.01740797 0.014858477 0.9809505
T3-T1 -0.009002020 -0.02513524 0.007131204 0.3864971
T3-T2 -0.007727273 -0.02457788 0.009123330 0.5252045
beta_anova(beta_div_T1_T2_T3, Herbicide = "Atrazine-Mesotrione")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.02773 0.013867 11.41 2.22e-05 ***
Residuals 173 0.21026 0.001215
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 -0.007369697 -0.02308574 0.008346347 0.5100816
T3-T1 -0.028906061 -0.04395303 -0.013859094 0.0000310
T3-T2 -0.021536364 -0.03658333 -0.006489397 0.0025367
beta_anova(beta_div_T1_T2_T3, Herbicide = "Glyphosate")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.02597 0.012985 14.9 1.33e-06 ***
Residuals 142 0.12374 0.000871
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 -0.02068822 -0.03474249 -0.006633939 0.0018742
T3-T1 -0.03156970 -0.04562397 -0.017515421 0.0000012
T3-T2 -0.01088148 -0.02562173 0.003858768 0.1909546
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide
index <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Handweeded", "Non-Treated")
sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])
#+ scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T1_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Mode")
T1_beta_chemical_combined_plot <- T1_beta_chemical_combined$plot
T1_beta_chemical_combined_plot<- T1_beta_chemical_combined_plot + theme_classic() + guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T1_beta_chemical_combined_plot
Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing missing values (geom_point).
my_legend <- get_legend(T1_beta_chemical_combined_plot)
Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing missing values (geom_point).
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_beta_combined_legend.pdf")
Saving 7.29 x 4.51 in image
T1_beta_chemical_combined_plot<-T1_beta_chemical_combined_plot+ theme(legend.position = "none")
T1_beta_chemical_combined_plot
Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing missing values (geom_point).
T2_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Mode")
T2_beta_chemical_combined_plot <- T2_beta_chemical_combined$plot
T2_beta_chemical_combined_plot<- T2_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T2_beta_chemical_combined_plot
T3_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Mode")
T3_beta_chemical_combined_plot <- T3_beta_chemical_combined$plot
T3_beta_chemical_combined_plot<- T3_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T3_beta_chemical_combined_plot
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
ggarrange(T1_beta_chemical_combined_plot, T2_beta_chemical_combined_plot, T3_beta_chemical_combined_plot, ncol = 1)
Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta_chemical_combined.pdf", width = 5, height = 10)
T1_beta_df_chemical_combined <- T1_beta_chemical_combined$data
T2_beta_df_chemical_combined<- T2_beta_chemical_combined$data
T3_beta_df_chemical_combined<- T3_beta_chemical_combined$data
T1_beta_df_chemical_combined$Time<-"T1"
T2_beta_df_chemical_combined$Time<-"T2"
T3_beta_df_chemical_combined$Time<-"T3"
m1<-aov(beta_div_value ~ group, data = T1_beta_df_chemical_combined)
summary(m1)
TukeyHSD(m1, "group")
boxplot(beta_div_value ~ group, data = T1_beta_df_chemical_combined)
m2<-aov(beta_div_value ~ group, data = T2_beta_df_chemical_combined)
summary(m2)
TukeyHSD(m2, "group")
boxplot(beta_div_value ~ group, data = T2_beta_df_chemical_combined)
m3<-aov(beta_div_value ~ group, data = T3_beta_df_chemical_combined)
summary(m3)
TukeyHSD(m3, "group")
boxplot(beta_div_value ~ group, data = T3_beta_df_chemical_combined)
beta_div_chemical_combined_T1_T2_T3 <- rbind(T1_beta_df_chemical_combined, T2_beta_df_chemical_combined, T3_beta_df_chemical_combined)
beta_anova(beta_div_chemical_combined_T1_T2_T3, Herbicide = "Chemical")
beta_anova(beta_div_chemical_combined_T1_T2_T3, Herbicide = "Hand")
beta_anova(beta_div_chemical_combined_T1_T2_T3, Herbicide = "Non-Treated")
treatment to control
plotDistances = function(p, m, s, d) {
# calc distances
wu = phyloseq::distance(p, m)
wu.m = melt(as.matrix(wu))
# remove self-comparisons
wu.m = wu.m %>%
filter(as.character(Var1) != as.character(Var2)) %>%
mutate_if(is.factor,as.character)
# get sample data (S4 error OK and expected)
sd = data.frame(sample_data(p)) %>%
select(s, d) %>%
mutate_if(is.factor,as.character)
sd$Herbicide <- factor(sd$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
# combined distances with sample data
colnames(sd) = c("Var1", "Type1")
wu.sd = left_join(wu.m, sd, by = "Var1")
colnames(sd) = c("Var2", "Type2")
wu.sd = left_join(wu.sd, sd, by = "Var2")
#remove this line to plot all comparisons.
#wu.sd = wu.sd %>% filter(Type1 == "Hand" | Type1 == "Non-Treated")
# plot
ggplot(wu.sd, aes(x = Type2, y = value)) +
theme_bw() +
geom_point() +
geom_boxplot(aes(color = ifelse(Type1 == Type2, "red", "black"))) +
scale_color_identity() +
facet_wrap(~ Type1, scales = "free_x") +
theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
ggtitle(paste0("Distance Metric = ", m))
}
a<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T1"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
a <- a + ggtitle("Time 1 Bray-Curtis Dissimlarities")
#ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_allgroup_beta.pdf")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_allgroup_beta_multicomparison.pdf")
b<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T2"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
b <-b + ggtitle("Time 2 Bray-Curtis Dissimlarities")
#ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_allgroup_beta.pdf")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_allgroup_beta_multicomparison.pdf")
c<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T3"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
c<- c + ggtitle("Time 3 Bray-Curtis Dissimlarities")
#ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_allgroup_beta.pdf")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_allgroup_beta_multicomparison.pdf")
library(ggpubr)
ggarrange(a, b, c, ncol = 1)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_allgroup_beta_multi_comparison.pdf", width = 7, height = 10)
Taxon abundance bar plot
#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
"#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
"#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
"#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
"#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
"#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
"#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
"#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
"#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
"#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
"#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
"#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
"#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
"#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
"#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
"#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
"#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
"#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
"#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
"#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
"#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
"#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
"#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
"#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",
"#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
"#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
"#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
"#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
"#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
"#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
"#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
"#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",
"#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
"#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")
phylumGlommed <- tax_glom(ps_rare, "Phylum")
#t1
phylumGlommed_herb_t1 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T1"), group = "Herbicide")
phylumGlommed_herb_t1 <- transform_sample_counts(phylumGlommed_herb_t1, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t1)$Herbicide <- factor(sample_data(phylumGlommed_herb_t1)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t1, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_Taxon_barplot_t1.pdf")
#t2
phylumGlommed_herb_t2 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T2"), group = "Herbicide")
phylumGlommed_herb_t2 <- transform_sample_counts(phylumGlommed_herb_t2, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t2)$Herbicide <- factor(sample_data(phylumGlommed_herb_t2)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t2, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/16S_Taxon_barplot_t2.pdf")
#t3
phylumGlommed_herb_t3 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T3"), group = "Herbicide")
phylumGlommed_herb_t3 <- transform_sample_counts(phylumGlommed_herb_t3, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t3)$Herbicide <- factor(sample_data(phylumGlommed_herb_t3)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t3, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/16S_Taxon_barplot_t3.pdf")
Combined herbicide and time bar plot for exploration
sample_data(ps_rare)$herb_time<-paste(sample_data(ps_rare)$Herbicide, sample_data(ps_rare)$Time, sep = "_")
ps_rare_for_barplot <- prune_taxa(taxa_sums(ps_rare) > 50, ps_rare)
plot_bar(ps_rare_for_barplot, x= "herb_time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_BarPlot_Herbicide_Time.pdf", width = 20, height = 11)
Linear modeling of abundant taxa.
Tax_glom_Subset <- function (physeq, y = "taxLevel", nreturns = "Number of returns"){
ps_1<- tax_glom(ps_rare_sub, taxrank = y )
myTaxa <- names(sort(taxa_sums(ps_1), decreasing = TRUE)[1:nreturns])
ps_1_sub <- prune_taxa(myTaxa, ps_1)
return(ps_1_sub)
}
ps_rare_family_top25<-Tax_glom_Subset(physeq = ps_rare, nreturns = 25, y = "Family")
myTaxa <- names(sort(taxa_sums(ps_rare), decreasing = TRUE)[1:25])
ps_rare_asv_top25 <- prune_taxa(myTaxa, ps_rare)
#explore top 25 taxa with plot bar
plot_bar(ps_rare_family_top25, x= "herb_time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
plot_bar(ps_rare_family_top25, x= "Time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
plot_bar(ps_rare_family_top25, x= "Herbicide", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
#write function to wrangle data prior to anova
abund_aov_wrangle <- function (physeq, y = "Tax_Level"){
tax<-tax_table(physeq)[,y]
meta<-data.frame(sample_data(physeq))
counts<-data.frame(otu_table(physeq))
rownames(counts) <- tax[,1]
counts<-data.frame(t(counts))
counts$Time <- meta$Time
counts$Herbicide <- meta$Herbicide
counts$Herb_time <- meta$herb_time
return(counts)
}
test<-abund_aov_wrangle(ps_rare_family_top25, y = "Family")
mod_abund<-function(count_tab, IV = "Groups to be tested") {
for(j in 1:length(unique(count_tab[,"Herbicide"]))){
data <- count_tab %>% filter(Herbicide == unique(count_tab$Herbicide)[j])
#change this to the number of returns from the tax_glom_subset function
for (i in 1:25) {
mod <- aov(unlist(data[i]) ~ matrix(data[,IV]))
#sanity check
#print(c(j, i))
if(summary(mod)[[1]][["Pr(>F)"]][1] <= 0.05) {
#print(summary(mod))
print(c(as.character(unique(count_tab[,"Herbicide"]))[j], names(data)[i]))
boxplot(unlist(data[i]) ~ unlist(data[IV]), main =paste(names(data[i]), as.character(unique(count_tab[,"Herbicide"]))[j]), xlab= "Time", ylab="Abundance")
}
}
}
}
mod_abund(test, IV = "Time")